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7 Abstract 

8 In a recent paper we presented a new ultra efficient numerical method for solving 

9 kinetic equations of the Boltzmann type |18) . The key idea, on which the method 

10 relies, is to solve the collision part on a grid and then to solve exactly the transport 

11 part by following the characteristics backward in time. On the contrary to classical 

12 semi-Lagrangian methods one does not need to reconstruct the distribution function 

13 at each time step. This allows to tremendously reduce the computational cost and to 

14 perform efficient numerical simulations of kinetic equations up to the six dimensional 

15 case without parallclization. However, the main drawback of the method developed 

16 was the loss of spatial accuracy close to the fluid limit. In the present work, we 

17 modify the scheme in such a way that it is able to preserve the high order spatial 

18 accuracy for extremely rarefied and fluid regimes. In particular, in the fluid limit, 

19 the method automatically degenerates into a high order method for the compressible 

20 Euler equations. Numerical examples are presented which validate the method, show 

21 the higher accuracy with respect to the previous approach and measure its efficiency 

22 with respect to well known schemes (Direct Simulation Monte Carlo, Finite Volume, 

23 MUSCL, WENO). 

24 Keywords: Kinetic equations, discrete velocity models, semi Lagrangian schemes, Boltzmann- 

25 BGK equation, Euler solver, high-order scheme. 



26 1 Introduction 

27 The kinetic equations provide a description of non equilibrium gases and more generally 

28 of particle systems [21 [9]. The distribution function, which describes the evolution of the 

29 system, depends, in the most general case, on seven independent variables: the time, the 
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30 physical and the velocity space. It turns out that the numerical simulation of these kind 

31 of equations with deterministic techniques presents several drawbacks due to the large 

32 dimension of the problem. On the other side of the spectrum of the numerical techniques 

33 used to approximate kinetic equations, there are the probabilistic methods HI [71 [3 132] , 

34 These methods and in particular Monte Carlo methods (DSMC) are extensively used due 

35 to their very low computational cost especially in the multidimensional cases, compared 

36 to finite volume, finite difference or spectral methods [211 ESI [30l [HI |33l |38l [Ml [37] , 

37 However, DSMC furnishes only poorly accurate and fluctuating solutions which cannot 

38 be easily ameliorated. This is especially true in non stationary situations in which time 

39 averages techniques turn to be useless. 

40 Many different works have been dedicated to reduce some of the disadvantages of 

41 Monte Carlo methods. We quote [7j for an overview on efficient and low variance Monte 

42 Carlo methods. Let us remind to the works of Homolle and Hadjiconstantinou [26] and 

43 [27] and of Dimarco, Pareschi and Degond [HI [20l [HI [IT] for some applications of variance 

44 reduction techniques to kinetic equations in transitional and general regimes. We recall 

45 also the works of Boyd and Burt [6] and of Pullin [39] who developed a low diffusion 

46 particle methods for simulating compressible inviscid flows. 

47 In this work, we continue the development of a new deterministic ultra fast method 

48 which permits to solve kinetic equations of the Boltzmann type [18]. The scheme is 

49 based on the classical discrete velocity models (DVM) approach [U [33l [Ml [37]. The 

50 DVM models are obtained by discretizing the velocity space into a set of flxed discrete 

51 velocities [H [Ml [Ml [M]- As a result of this discretization, the original kinetic equation 

52 is then represented as a set of linear transport equations plus a source term. The source 

53 term describes the collisions or the interactions between the particles and couples all the 

54 equations of the resulting system. In order to solve the transport part of the DVM model, 

55 many different techniques can be employed like flnite difference, flnite volume or fast 

56 methods [H [22l [211 [311 [Ml [25] . One of the most common strategies for solving this kind of 

57 problems is the semi Lagrangian approach [10\ \TT\ [22l \^\ which will be also the basis of 

58 the method here developed. Unfortunately, for each of the method cited, we recall that, 

59 the computational effort needed for solving the full six dimensional equation, prevents still 

60 nowadays realistic simulations even with parallel machines. 

61 To overcome the problem of the excessive computational cost, we recently proposed 

62 in [l8j to use a splitting method to separate the transport from the collision step [131 142j . 

63 Then, we used a Lagrangian technique which exactly solves the transport step on the entire 

64 domain and we projected the solution on a grid to compute the contribution of the collision 

65 operator. The resulting scheme (Fast Kinetic Scheme, FKS) shares many analogies with 

66 semi-Lagrangian methods \10\ [TT\ \T2\ [22l [JO] and with Monte Carlo schemes [29], but 

67 on the contrary to them, the method is as fast as a particle method while the numerical 

68 solution remains fully deterministic, which means that there is no source of statistical 

69 error. When used to solve the limiting fluid model, the FKS method shares also some 

70 analogies with the so-called Lattice Boltzmann methods [1], but on the contrary to them 

71 its application is not limited to dense flows, all the regimes from rarefied to dense can be 

72 studied with such approach. Thanks to this new scheme, we were able to compute the 

73 solution of the full six dimensional kinetic equation on a laptop for acceptable mesh sizes 
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74 and in a reasonable amount of time (about ten hours for 100'^ space x 12^ velocity space 

75 mesh points for 110 time steps up to i = 0.1 on the spherical Sod- like problem). 

76 However, the method developed, exhibited some limitations in term of spatial accuracy. 

77 In fact, it was only poorly accurate when used to compute solutions close to the fluid limit. 

78 We showed in [TH] that, in these cases, the computed solutions laid between a first order 

79 and a second order MUSCL scheme. In the present work, we developed a strategy which 

80 permits to preserve a desired high-order of accuracy in space for all the different regimes 

81 from the extremely rarefied to the high dense cases. In particular, it permits to recover 

82 the solution of the compressible Euler equations, when the number of collisions tends to 

83 infinity, with an high order shock capturing scheme. The modification introduced consists 

84 in coupling the fast kinetic scheme (FKS) to a solver for the compressible Euler equations, 

35 then to match the moments obtained from the solution of the macroscopic equations with 

36 those obtained from the solution of the equilibrium part of the kinetic equation. Finally, 

87 the solution, in term of the moments, is recovered as a convex combination of the two 

88 contributions: the macroscopic and the microscopic parts. We will show that the in- 
39 troduction of a macroscopic solver will not increase dramatically the computational cost, 

90 instead this modification will represent only a fraction of the time employed for computing 

91 the solution. In this work, the interaction term between the particles is the BGK collision 

92 operator |23j. However, the high order version of the Fast Kinetic Scheme can in principle 

93 be extended to other collisional operators as the Boltzmann one. 

94 

95 The article is organized as follows. In section [21 we introduce the Boltzmann-BGK 

96 equation and its properties. In section [3l we present the discrete velocity model (DVM). 

97 Then, in section U we present the fast kinetic scheme (FKS) and the High Order Fast 
93 Kinetic Scheme (HOFKS). Several test problems which demonstrate the accuracy and the 

99 strong efficiency of the new method are presented and discussed in section \E\ Some final 

100 considerations and future developments are finally drawn in the last section. 

101 2 Boltzmann-BGK Equation 

102 The equation to be solved is the following: 

dtf + vV,f = ^iMf-f), (1) 

103 with the initial condition /(x, v,t = 0) = fo{x, v). The non negative function / = /(x, v, t) 

104 describes the time evolution of the distribution of particles which move with velocity 

105 V € TZ'^ in the space x € O C TZ'^ at time t > 0. For simplicity, in the description of the 

106 method we will do the hypothesis that the dimension of the physical space is the same of 

107 the dimension of the velocity space d. However, the method is not restricted to this par- 

108 ticular choice and it is possible to consider different dimensions between the space and the 

109 velocity in order to obtain different simplified models. The type of interactions term which 
no characterizes the kinetic equation in ([1]) is the so-called BGK relaxation operator. With 
111 this choice the collisions are modeled by a relaxation towards the local thermodynamical 
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112 equilibrium defined by the Maxwellian distribution function M f 



M, = M,[p,.,r](.) = (^exp [-^^ j , (2) 

113 where p ^ IZ, p > and u G TZ'^ are respectively the density and mean velocity while 

114 6 = RT with T the temperature of the gas and R the gas constant. The macroscopic 

115 values p,u and T are related to / by: 

fdv, u = - I vfdv, ^ = ~:i I b ~ u\'^fdv, (3) 
T^d p J-jid pd Jfid 

116 while the energy E is defined by 

E=^I^Jvffdv = ^p\u\^ + ^pe, (4) 

117 The parameter r > in ([1]) is the relaxation time. We refer to section [5] for the numerical 

118 values chosen. 

119 Formally, when the number of collision goes to infinity, which means r — >■ 0, the 

120 function / converges towards the Maxwellian distribution. In this limit, if we consider the 

121 BGK equation ([T]) and we multiply it by 1, v, and then we integrate with respect 

122 to f , we get the so-called Euler system of compressible gas dynamics equations 

^ + V,-{pu) = 0, 
dpu „ 

— + V.,-i{E+p)u) = 0, 
p = p9, E ='^p9 + ]^p\u\^. 

123 In the following, we will combine the BGK equation ([1]) with the compressible Euler 

124 equation ^ to get our High Order Fast Kinetic scheme. 



125 3 The Discrete Velocity Models (DVM) 

126 The principle of Discrete Velocity Model (DVM) is to set a grid in the velocity space and 

127 thus to transform the kinetic equation ([T]) in a set of linear hyperbolic equations with 

128 source terms. We refer to the works of Platkowski |36] and of Mieussens [30] for the 

129 description of this approach and we remind to them for the details. In the following, we 

130 will briefly describe the idea and we will introduce the notations which will be used. 

131 Let /C be a set of M multi-indices of N'^, defined hy K, = {k = {k^^)f^^, feW < iCW}, 

132 where {K^^^} are some given bounds. We introduce a Cartesian grid V of M*^ by 

V = {ufc = kAv + a, keIC}, (6) 
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133 where a is an arbitrary vector of M and Aw is a scalar which represents the grid step in 

134 the velocity space. We denote the discrete collision invariants on V by = (1, u^, 

135 Now, in this setting, the continuous distribution function / is replaced by a A^— vector 

136 f)c{x,t), where each component is assumed to be an approximation of the distribution 

137 function / at location t;^: 

fK{x,t) = {fk{x,t))k, fk{x,t) ^ f{x,Vk,t). (7) 

138 The fluid quantities are then obtained from thanks to discrete summations on V: 

U{x,t) = ^mkfkix,t) Av = {mkfk{x,t))x.- (8) 
k 

139 The discrete velocity BGK model consists of a set of A'^ evolution equations for of the 

140 form 

dtfk + Vk ■ y.fk = -{£k[U] - fk), k = l, .., N (9) 
r 

141 where ffc[f7] is a suitable approximation of defined next. The DVM approach deserves 

142 some remarks. 

143 Remark 1 



144 • When dealing with discrete velocity methods, one needs to truncate the velocity space 

145 and to fix some bounds. This gives the number N of evolution equations Of 

146 course, the number N is chosen as a compromise between the desired precision in 

147 the discretization of the velocity space and the computational cost, while the bounds 

148 are chosen to give a correct representation of the flow. This implies that the dis- 

149 Crete velocity set must be large enough to take into account large variations of the 

150 macroscopic quantities which may appear during the evolution of the problem. On 

151 the other hand, the number of mesh points should be sufficiently large to guarantee 

152 that the small variations of the macroscopic quantities are well described. 

153 • The exact conservation of macroscopic quantities is impossible, because in general 

154 the support of the distribution function is non compact. This is the case for instance 

155 of the Maxwellian equilibrium distribution. Thus, in order to conserve macroscopic 

156 variables, different strategies can be adopted, two possibilities are described in 123[ 

157 \^(J^ . Moreover, the approximation of the equilibrium distribution Mf by £k\U] must 

158 be carefully chosen in order to satisfy conservations of physical quantities. In the 

159 following section we will discuss our choices in details. 



160 4 Fast Kinetic Schemes (FKS) and High Order Fast Kinetic 
Schemes (HOFKS) 

162 In this section we recall the FKS method and then we will introduce a new class of schemes 

163 which enables to get high order spatial accuracy (HOFKS). Before, we will discuss and 

164 propose a solution for the problem of lack of conservation of the macroscopic quantities 
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165 which characterizes the class of discrete velocity models we are dealing with. We first 

166 of all introduce a Cartesian uniform grid in the physical space, FKS schemes are in fact 

167 based on uniform meshes. The extension of this class of methods to general meshes is not 

168 trivial but nevertheless under study. The mesh is defined by the set J oi N multi-indices 

169 of N'^, which is ^ = {j = J^*^ < J^^^}-, where {J^*^} are some given bounds which 

170 represent the boundary points in the physical space. The grid X of is then given by 

X = {xj = jAx + b, j G J}, (10) 

171 where d represents at the same time the dimension of the physical space and the dimension 

172 of the velocity space. The form of the physical space is determined by the vector b of 

173 and Ax is a scalar which represents the grid step in the physical space. We consider a 

174 third discretization which is the time discretization = nAt with n S N. We will at the 

175 end of the section discuss the time step limitations and the CFL condition. 

176 4.1 Conservative Discrete Velocity Models (DVM) 

177 Suppose a continuous in phase space distribution function is given, i.e. f{xj,v,t"'), with 

178 moments U{xj,t") for every j G J and n > 0. We proceed into two steps. First we define 

f^^j, = f{xj,Vk,tn), (11) 

179 which is the pointwise distribution value in phase space. Observe that, due to the trunca- 

180 tion of the velocity space and to the finite number of points with which / is discretized, the 

181 moments of /j^^ differ from the original moments U{xj,t"'). In fact the discrete moments 

182 of this distribution are 

U^ = {mkjf,,Av)jc^U^, p;;<P^, 0^<0]. (12) 

183 Different strategies can be adopted to restore the correct moments. Our choice, which is 

184 the second step of the conservative DVM model, consists in defining the approximated 

185 distribution function as the distribution closer in the discrete L2 norm to = 

186 f{xj,Vk,tn) for which the moments are exactly the macroscopic quantities we want to 

187 preserve , i.e 

Uj = {mUlu^v)^. (13) 

188 In order to find this distribution we make use of a simple constrained Lagrange multi- 

189 plier method |23] . where the constraints are mass, momentum and energy of the solution. 

190 Let us recall the technique from (23]. For each spatial cell, let define the pointwise distri- 

191 bution vector 

^" ••'t^)^ (14) 

192 let also define the vector containing the corrected distribution which fulfills the conserva- 

193 tion of moments we are searching for 
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194 and the matrix which contains the discretization parameters C G ]r(^+2)x^_ At this point, 

195 conservation can be imposed in each cell and at any time index n solving the following 

196 constrained optimization problem: 

Given S M^, C G R(.d+2)xN ^ jjn ^ ^{d+2)xl^ 

find f/ G such that (16) 
II/" — /" II 2 is minimized subject to the constraint C/" = C/". 

197 Thus, let A G W^'^'^ be the Lagrange multiplier vector, the objective function to be 

198 minimized, in each cell, is given by 

N 

m^^) = E i& - + ^^(^/r - f^D- (17) 

k=l 

199 The above equation can be solved explicitly. The searched distribution function is then 

200 is ^ 

f^ = f^ + c^{cc^rHu^-cf^), yj^J. (18) 

201 We end this part defining the approximated equilibrium distribution [[/"] , or equiva- 

202 lently S^f^[U]. The discretization of the Maxwellian distribution Mf{x, v,t), should satisfy 

203 the same properties of conservation of the distribution /, i.e. 

[/; = (mfc/;;, Av),c = [rnkSkP^] Av),c. (19) 

204 To this aim, observe that the natural approximation £"/;[[/"] = Mf{xj,Vk,tn) = Mj[C/"] 

205 cannot satisfy these requirements. Thus, the calculation carried out above for the defi- 

206 nition of the approximated distribution /, should also be performed for the equilibrium 

207 distribution Mj. This should be done each time we invoke the equilibrium distribution 

208 during the computation as explained in the next subsection. The function £\U] is therefore 

209 given by the solution of the same minimization problem defined in ()16p . and its explicit 

210 value is given mimicking (llSp by 

E[UJ]= Mf[U]'] + C^{CC'^)-\Uj -CMf\Ul^]), "ij^J. (20) 

211 Notice that the computation of the new distributions / and E only involves a matrix-vector 

212 multiplication. In fact, matrix C only depends on the parameter of the discretization and 

213 thus it is constant in time. In other words matrices C and C"^(CC"^)~^ can be precomputed 

214 and stored in memory while initializing the problem. 

215 Remark 2 

216 • For FKS schemes, we need to solve the above minimization problem for the initial 

217 data f{xj,Vk,t^ = 0) and for the distribution £[Uj'] at each time index n. In fact, 

218 once the conservation is guaranteed for f for t = 0, this is also guaranteed for the 

219 entire computation because the exact solution is used for solving the transport step. 



7 



220 • The only possible source of loss of conservation for this type of schemes is due to 

221 the solution of collision term and to the way in which the equilibrium distribution is 

222 discretized. This aspects will be detailed in the next subsection. 

223 • The conservation technique described in il6\} does not assure the positivity of the 

224 distribution function /"^. It may happen in fact that during the constrained min- 

225 imization procedure becomes negative for some values of k. In practice, we did 

226 not observe this phenomenon to create instability in the solution. However, in the 

227 cases in which positivity is strictly demanded, as for instance for the full Boltzmann 

228 operator discretized with spectral methods, alternative techniques should be designed. 



229 4.2 FKS schemes 

230 The main features of the method developed in [18] can be summarized as follows: 

231 • The BGK equation is discretized in velocity space by using the DVM model. The 

232 distribution / as well as the Maxwellian Mf are initialized by the conservative DVM 

233 method detailed in Section 14. 1[ 

234 • A time splitting procedure is employed between the transport and the relaxation op- 

235 erators for each of the resulting N evolution equations ([9]) . First order time splitting 

236 is considered. In principle, others more sophisticated splitting can be employed [42] . 

237 • The transport part is solved exactly and continuously in space, this means that no 

238 spatial mesh is involved. The initial data of this step is given by the solution of the 

239 relaxation operator. 

240 • The relaxation part is solved on the spacial grid. The initial data for this step is given 

241 by the value of the distribution function in the center of the cells after the previous 

242 transport step. Each time the equilibrium distribution is invoked, conservation is 

243 retrieved through equation (I20p . We need to impose conservation of the macroscopic 

244 quantities for the equilibrium distribution only. 

245 Let us give the details of the method. We recall that the FKS methods are constructed 

246 on uniform grids. Let /^^ be the pointwise initial data, solution of equation (I18p with 

247 f^^. = f{xj,Vk,t = 0). Let also £jj^[U] be the initial equilibrium distribution solution 

248 of equation (pO|) with M?^ = Mf{xj,Vk,t = 0) defined at points Xj at i = as the 

249 distribution /. We start describing the first step of the method starting at t*^ = 0. 

250 The scheme is then generalizable to the generic time step [t"'; i*^"*"^] starting from t^. 



251 First time step Let us describe the transport and relaxation stages. 

252 Transport stage. We solve N linear transport equations of the form: 

dtfk + Vk-V,fk = 0, k = l,...,N. (21) 
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253 The idea of FKS schemes is to solve the transport part continuously in space instead 

254 of solving it only on the mesh points. To this aim, we define for each of the N 

255 equations a piecewise constant function in space as 

7fc(a;,t0 = 0) = /0fc VxG [x,-„i/2,x,+i/2], k = l,...,N. (22) 

256 Now, the exact solution of the N equations at time = + At = At is given by 

Tk{x)=J{x-VkAt), k = l,...,N. (23) 

257 Observe that, with this choice, we do not need to reconstruct the distribution func- 

258 tion / as for instance in the semi-Lagrangian schemes [21^ 122): the shape of the func- 

259 tion in space is in fact known at the beginning of the computation and it remains 

260 so through the duration of the computation. To be more precise the distribution 

261 function is transported in time with constant velocity so the discontinuities remain 

262 at the same relative locations. It remains also a piecewise constant function. The 

263 relaxation step, as finite difference methods, is solved only on the grid points. This 

264 means that only the value of the distribution function / and the macroscopic quan- 

265 titles in the centers of the cells are needed for this step. From the exact solution of 

266 the function we can immediately recover these values at the cost of one simple 

267 vector multiplication. 

268 Relaxation stage. This step is local to the grid, this means that we solve the following 

269 ordinary differential equation: 

dtfj,k = ^{£jAU]-fj,k), k = l,...,N, j = l,...,M, (24) 

270 where the initial datum is the result of the transport step at points xj at time 

271 t^ =t° + At 

Tk{xj) = Rxj-VkAt), k = l,...,N, j = l,...,M. (25) 

272 To solve equation (j25p we need the value of the equilibrium distribution £ at the 

273 center of the cell after the transport stage. In order to compute the Maxwellian, 

274 the macroscopic quantities in the center of the cells, i.e. the density, the mean 

275 velocity and the temperature, are given by summing the local value of the discrete 

276 distribution / over the velocity set: {mkf* f^Av)ic = U*, for all j = 1, . . . , M, where 

277 /*^ = f*k{xj). The discrete equilibrium distribution at time t^ = t^ + At, £* = Sji^, 

278 is the solution of equation (|2U|) with moments U* = Uj, for all j = 1, . . . , M. Observe 

279 in fact that, the Maxwellian distribution does not change during the relaxation step. 

280 In other words during this step the macroscopic quantities remain constants. We 

281 can now compute the solution of the relaxation stage as 

fl, = eM-^t/T)fl, + (1 - eM-^t/T))£l,[U]. (26) 

282 The above equation furnishes only the new value of the distribution / at time t^ = 

283 t^ + At = At in the center of each spatial cell for each velocity • However in order 
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284 to continue the computation, we need the value of the distribution / in all points of 

285 the space. Let us assume that the equilibrium distribution Mf has the same shape 

286 than the distribution / in space. Thus, starting from the pointwise value of £ we 

287 define a piecewise constant function in space £k for each velocity Vk- The values of 

288 this piecewise constant function are the values computed in the center of the spatial 

289 cells. In other words, one defines 

£*kix) = £k{x,t^) = Slk, Vx suchthat /^(j;) = /^(rEj), j = 1, . . . , M. (27) 

290 We can now rewrite the relaxation term directly in term of spatial continuous func- 

291 tion as 

Jlix) = J,{x, At) = exp(-At/T)7*(x) + (1 - exp(-At/r))^*(x)[C/]. (28) 

292 For each velocity Vk the original shape in space for the distribution fk is preserved 

293 throughout the computation, and, as a consequence it drastically reduces the com- 

294 putational cost because no reconstruction is needed. 

295 The time marching procedure can be now be described. 

296 Generic time step Given the value of the distribution function /^(x), for 

297 all k = 1, . . . ,N, and all x G M'^ at time t^, the value of the distribution at time t^~^^, 

298 (x), is given by 

Tk(x)=Tk{x-VkAt), k = l,...,N (29) 

299 

7^+'(x) = exp(-At/r)7;(x) + (1 - exp(-Ai/T))^r'(^)[^], A; = 1, . . . ,iV, (30) 

300 where f)!^^(x)[C/] is a piecewise constant function with the discontinuities located in the 

301 same positions as the distribution It is computed considering the solution of the 

302 minimization problem (|20p relative to the moments value in the center of each spatial cell 

303 after the transport stage: U^~^^, j = 1, • • • ,M. These moments are given by computing 

304 {mj^f^ f^Av)]c where /*^ is the value that the distribution function takes after the transport 

305 stage in the center of each spatial cell. 

306 Remark 3 



307 • Due to the fact that the relaxation stage preserves the macroscopic quantities, the 

308 scheme is globally conservative by construction. In fact, at each time step, the change 

309 of density, momentum and energy is only due to the transport step. This latter being 

310 exact, does preserve the macroscopic quantities as well as the distribution function. 

311 • For the same reason, except for the constrained optimization procedur^, the scheme 

312 is also unconditionally positive. More precisely, if f]^ {x) > 0, and k = 1, . . . , M and 



Observe that the positivity of the constrained optimization step can be forced introducing an inequality 
constraint of the type f^^. > or S"^. > 0. However, the introduction of such a step will cause the 
minimization step to be solved numerically instead of analytically. This will means that the computational 
cost of the method will increase. In the present work, we did not attack this problem and we remind to 
the future for the development of strictly positivity preserving fast kinetic schemes. 
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313 the optimization procedure preserves positivity, then {x) > if the initial datum 

314 is positive f^{x) > for all k = 1, . . . , M. In fact, the transport maintains the shape 

315 of f unchanged in space while the relaxation towards the Maxwellian distribution is 

316 a convex combination of Mf and f{x — VkAt) both being positive. 

317 • The time step At is constrained by the CFL condition 

At max ( ) < 1 = CFL. (31) 

318 Observe that this choice is not mandatory, in fact the scheme is stable for every 

319 choice of the time step. However, being based on a time splitting technique the error 

320 is of the order of At. This suggests to take the usual CFL condition in order to 

321 maintain the time splitting error small enough. 

322 • Some experiments have been done on the influence of the CFL condition on the 

323 schemes. The results showed that, for the cases tested, up to CFL = 5 in \31\) . the 

324 FKS scheme provides a solution very close to the case CFL = 1 for all the values 

325 of the Knudsen number. Moreover, when the Knudsen number is large, i.e the BCK 

326 equations are very close to a free transport equation, using larger values of the CFL 

327 number does not cause any more degradation to the global accuracy because de facto 

328 the FKS scheme solves the transport term exactly. 



329 4.3 HOFKS schemes 

330 As observed in [18] the FKS scheme performs very well in collisionless or almost collisionless 

331 regimes. In these cases, in fact, the relaxation stage is neglectable and only the exact 

332 transport does play a role. However, when moving from rarefied to dense regimes the 

333 projection over the equilibrium distribution becomes more important. Thus, the accuracy 

334 of the scheme was expected to diminish in fluid regimes, because the projection method 

335 is only first order accurate. These behaviors were, in fact, observed in the numerical 

336 simulations performed [18]. In the present paper, we developed a method which preserves 

337 the high spatial accuracy observed with the FKS schemes for rarefied regimes and which 

338 becomes a high order shock capturing scheme applied to the kinetic equation ([TJ in the fluid 

339 limit. This means that throughout all possible regimes, from fluid to extremely rarefied 

340 flows, the new scheme maintains high accuracy in space. Moreover, the new method does 

341 not cause the computational efficiency to drop down. As shown in the numerical test 

342 section the high order fast kinetic schemes (HOFKS) still works with computational costs 

343 close to the original FKS method and, for unsteady problems, in which time averaging are 

344 unusable, it is still much faster than DSMC methods. We recall that we only focus on the 

345 spatial accuracy in this paper, we postpone to the future the development of high order 

346 schemes both in time and space. 

347 4.3.1 The general methodology 

348 The idea, onto which the method is based, is to solve the equilibrium part of the distribu- 

349 tion function with a macroscopic scheme instead of a kinetic scheme. In fact, observe that 
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350 at each time step, the relaxation stage consists in computing the distribution function /" 

351 as a convex combination of the transported distribution /* and a Maxwelhan distribution 

352 The Maxwellian distribution 8"^ is computed through the moments of /*. Then, in 

353 order to complete one step in time, the scheme solves the transport part which leads to 

354 the new intermediate distribution /* at time n + So now, in HOFKS scheme, we replace 

355 the moments at time index n + 1, obtained from the solution of the transport stage at time 

356 n with another set of moments. This new set of moments are computed through the same 

357 convex combination of the relaxation stage ()30p . However now, to the contrary of ()30p . 

358 the convex combination is performed between the moments which come from the solution 

359 of the transport part of the kinetic equation and the moments which are solutions of the 

360 compressible Euler equations. At this point, if the compressible Euler equations are solved 

361 with a high order shock capturing scheme in the limit r ^ the HOFKS corresponds to 

362 the same method for the macroscopic equations. We detail this new scheme in the sequel. 



363 In order to keep notations simple and compact we introduce three operators: the 

364 projection operator V, the relaxation operator T^At and the transport operator TXt which 

365 act on a time step At. 

366 From the kinetic variable / (or Mj) the projection operator computes the macroscopic 

367 averages U{xj,t'^) = U^, thus 

V]{f) = {m,f^^,Av)^ = {mA[U^]Av)K=Vj{£) = U{xj,n, (32) 

368 since the local Maxwellian Mj has the same moments of the distribution function /. 

369 The relaxation and transport operators solve the relaxation and transport steps for piece- 

370 wise constant functions and £k or equivalently for pointwise functions fj^k and 

371 £j^k[U] for all velocities Vk, k = 1, ..,N. The relaxation operator has the form 

7^A^(7) =A7+(1-A)^, (33) 

372 where A = exp(— At/r), whereas the transport operator reads 



TAt{f)=f{x-vAt). (34) 

373 In order to describe the HOFKS scheme, let us start, to the contrary of FKS scheme, from 

374 the relaxation step. Recall that starting either from the relaxation or from the transport 

375 step gives consistent splitting discretizations. Suppose that the distribution function (x) 

376 is known as well as the function f ^ (x) for all A; = 1, . . . , A^, then, as first step of the splitting 

377 we have 

r = T^At (T) = XT + (1 - A)r . (35) 

378 The distribution function /* is given by a convex combination of the transported distri- 

379 bution and the Maxwellian distribution. Then, the transport step, applied to the solution 

380 of the relaxation step, produces the so called kinetic solution (K) at time index n + 1 

17' = TAt (r) = TAt {xT) + TAt ((1 - x)r) . 
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381 On the other hand, the kinetic solution in terms of the macroscopic moments furnishes 

382 the following values 

C/x(x,,r+i) = rj+'(TAt{T)) (36) 

= vj+\TAt{xT)) + r;^^\TAt{{i-x)'r)) (37) 
= u,ixj,e+') + UM{xj,e+^). (38) 

383 In order to construct the HOFKS scheme we replace the moments UM{xj,t'^~^^) by the 

384 moments obtained solving the compressible Euler equations that we call UE{xj,t"~^^). 

385 (The details of the numerical scheme used will be given next.) 

386 Thus the final moments used in the solution are given by 

UHixj,e+') = C/,(x„r+i) + UE{xj,e+') (39) 

337 where Uh stands for hybrid. Before describing the last step which ensures consistency 

338 between the kinetic solver for the Maxwellian distribution and the macroscopic solver for 

389 the compressible Euler equations, we state the following result (see [191 [35] for a proof): 

390 Theorem 1 If we denote by UE{x,t + At) the solution of the Euler equations and 

391 with C/jvf(x,t + At) the solution of the kinetic equation in which we consider initial ther- 

392 modynamical equilibrium, i.e. f = S\U]. If in addition we consider as initial data 

393 UE{x,t) = UM{x,t) then 

Ue{x, t + At) = Um{x, t + At) + 0{At^). (40) 

394 

395 By virtue of the above result, we can replace the moments after the transport UM{xj.,t"') 

396 with UEixjjf^) at each time step without affecting the overall first order accuracy in time 

397 of the splitting method. However, to have consistency between the macroscopic solution 

398 and the kinetic discretized solution, it is necessary that the advected equilibrium satisfies 

T^r' {TAtiil - A)r)) = UE{x„t^+'), (41) 

399 namely the kinetic solution to the fluid equations in one time step should match the direct 

400 solution to the limiting fluid equations. This is not true in general. To solve this problem, 

401 we apply again the minimization method ()20p to find the new distribution Tai ^(1 — 

402 which shares the same moments than the macroscopic solution C/ij(xj,t"+^). Thus, we 

403 search for a distribution ^(1 — A)£'"^ which satisfies the following minimization prob- 

404 lem: 

Given rAt((l - X)£^) G M^, Ce M('^+2)><^, and [/^^(xj, t"+i) G m('^+2)x1, 

find r^((l - \)£f) G such that (42) 
||7At((l — A)£'") — 7^(((1 — A)£'")|[2 is minimized subject to the constraint 

c(TU{^-X)£^))=UE{xj,t^+'). 
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405 Then again, starting from the pointwise solution, we define as in the FKS schemes a new 

406 piecewise constant equihbrium function 7^4((1 — X)£){x) sharing its shape with / as 

TAt ((1 - A)?:) (x) = rAa(l - A)^:]:^) , Vxs.t. K(rE)=7:(x,), j = l,...,M. (43) 

407 Finally, the new distribution /, at time index n + 1, is defined as 

T''' = TAtiXT) + TUil - A)^), (44) 

408 while the new moments are given by (j39p . This somehow ends the methodology to design 

409 HOFKS schemes. What remains to be detailed is how the solution of the compressible 

410 Euler equations UE{xj,t"'~^^) = U^'^^ is computed. Remark that at this point any solver for 

411 the compressible Euler equations can be used. One example is proposed in the following. 



412 4.3.2 One example: MUSCL Finite Volume (FV) scheme 

413 As a first example we propose the MUSCL Finite Volume (FV) scheme. This is also the 

414 scheme used in the numerical test section. It reads, starting from U^j, 

^g-^g. , V'.+l/2(^g)-V'.-l/2(^g) „ 

At ^ Ax ' ^ ' 

415 where discrete fluxes are defined as in [28] by 

^Hi/2m) = limh) + nuij+i)) - ^c^m.+i - ui^) + ^(^j'^ - a;-) (46) 

416 with F{U) the flux of the compressible Euler equations and 

a]'^ = (F(f/^,,+i) ± aU^j+, - FiU^^^) ^ aU^^) ^{x]'^) (47) 

417 with if being the slope limiter, as instance we use the Van Leer slope limiter 

M = (48) 
1 + X 

418 where the variable is defined as following 



F{Ul._^,)±aUl._^,-F{Ul.)TaUi: ^ ' 



419 The above ratio of vectors is defined componentwise and a represent the eigenvalues of 

420 the Euler system. 
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421 4.3.3 Extensions and stability constraints 



We have proposed second order schemes to solve the compressible Euler equations. How- 
ever, any other solver can in principle be used, which could additionally increase the spatial 
accuracy of the HOFKS method, for instance WENO methods [H] or a genuine MOOD 
scheme [151 El • 

Finally, the time step At is chosen such that it satisfies the stability condition of the Euler 
solver, in fact we recall that the FKS is stable for all choices of the time step. This means 
that the time step is driven for the MUSCL scheme by 

At = -( . (50) 

with amax the largest eigenvalue of the Euler system. If a more CFL restrictive Euler 
solver has to be used, as instance with a WENO scheme, then the time step will be ac- 
cordingly reduced. 



433 5 Numerical tests 

434 In this section, we present several numerical tests to illustrate the main features of the 

435 method and the improvements with respect to the FKS scheme. The following methodol- 

436 ogy is adopted 



437 First, we test the HOFKS method on the one dimensional Sod shock tube. In this case, we 

438 compare two kinetic schemes (the FKS, a first order upwind scheme and the HOFKS 

439 a second order scheme) versus a finite volume (FV) upwind scheme, a second order 

440 finite volume MUSCL method (MUSCL) and a third order WENO finite difference 

441 scheme (WENO). For any fiuid regime the goal is threefold: (i) the two kinetic 

442 schemes produce valide results, (ii) the HOFKS is more accurate than FKS and, (iii) 

443 the accuracy of FKS lays in between FV and MUSCL and the accuracy of HOFKS 

444 lays in between MUSCL and WENO. 

445 In a second test case we use the exact smooth solution of the advected isentropic vortex 

446 of the 2D Euler equations to assess the effective numerical accuracy and rate of 

447 convergence of FKS, HOFKS and also the unlimited version of HOFKS. The solution 

448 being smooth an unlimited scheme can be used to measure the maximal accuracy 

449 that can be obtained with our choice of Euler solver. 

450 Then, in a third series of tests we solve a two dimensional-two dimensional BGK equation 

451 and we compare our method with a Monte Carlo scheme (DSMC) for r = 10~^, and, 

452 in the fluid limit, r = 10""^, with FV and MUSCL schemes. The goal is to show 

453 that a genuine multi-dimensional solutions with shocks and interaction of waves 

454 can be accurately captured with HOFKS in any fiuid regime. We also report the 

455 computational times for the two dimensional simulations for the HOFKS, the FKS 

456 and the DSMC methods. All simulations are performed on a mono-processor laptop 

457 machine. 
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458 5.1 ID Sod shock tube problem 

459 We consider the ID/ ID Sod test with 300 mesh points in physical and 100 points in 

460 velocity spaces. The boundaries in velocity space are set to —15 and 15. The left and 

461 right states are given by a density pL = I, mean velocity ul = and temperature = 5 

462 if < < 0.5, while pR = 0.125, ur = 0, Tr = A if 0.5 < x < 1. The gas is, at the 

463 initial state, in thermodynamical equilibrium. We repeat the same test with four different 

464 values of the Knudsen number, i.e. t = 10~^, r = 10~^, r = 5 x 10~^ and r = 10~^. 

465 We plot the results for the final time tgnai = 0.05 for the density (Figured]), the mean 

466 velocity (Figure [2]) and the temperature (figure [3]). In each figure we compare the HOFKS 

467 method with the FKS method. We reported also the solutions computed with a third order 

468 WENO method, a second-order MUSCL method and a first-order upwind method [28]. 

469 These numerical methods, used as reference, employ the same discretization parameters. 
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470 except for the time steps which for each scheme is chosen in order to satisfy the stabihty 

471 conditions. 

472 From Figures ([I]) to © we observe that the HOFKS, the FKS and the WENO methods 

473 give identical or almost identical results for r = 10~^, this result was expected. We build 

474 up the method in such a way that for larger r it behaves like the original fast kinetic 

475 scheme, because we knew that in these regimes the FKS already gave very good results. 

476 For larger values of r we found the same behaviors as for the case r = 10-2, thus we 

477 did not report simulations results. Starting from r = 10"^, some small differences arise 

478 between the HOFKS and the FKS methods, however both schemes are still very close to 

479 the WENO solution. For r = 5 x 10"^, we clearly see differences between the high order 
4B0 fast kinetic scheme and the fast kinetic scheme. In particular, we see that the HOFKS 

481 remains stick to the third order WENO scheme while the FKS not. This aspect is made 

482 very clear for r = 10"^. In this latter case, the HOFKS gives very good results while 
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4B3 the FKS scheme lays between a first and a second order space accurate method. The key 

484 point of the schemes developed is their very low CPU time consumption in comparison to 

485 other existing methods. This gain as expected is not so relevant for the one dimensional 
4B6 case, while it becomes very important for the two and the three dimensional cases. Thus, 
487 in the next subsection we report some two dimensional simulations together with their 
4B8 computational costs. 

489 5.2 2D isentropic vortex 

490 This vortex test case is a classical test to assess the accuracy of numerical schemes because 

491 this problem produces a genuine 2D smooth solution of Euler equations. The overall error 

492 produced by our kinetic schemes associates spacial discretization error, time discretization 

493 error and velocity discretization error. In our case the first order time discretization dooms 
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the overall scheme to remain first order accurate only. Nevertheless using a time step At of 
the order Ax^ or even smaller allows to reduce the time discretization error to a negligible 
quantity by respect to space discretization error. Doing so we can effectively measure the 
spacial accuracy of the overall scheme when the velocity error is kept small enough. 
To do so an isentropic vortex is introduced to a uniform mean flow, by small perturbations 
of velocity, density and temperature variables and is detailed in [5T], [33] as instance. The 
simulation domain Q is the square [0, 10] x [0, 10] and we consider an initial gas flow given 
by the following background condition poo = 1-0, u^o = 1-0, Voo = 1.0, Poo = 1-0, with 
a normalized ambient temperature = 1.0 computed with the perfect gas equation of 
state and 7 = 5/3. 




Exact density - 
Exact temperatjre 
Exact pressjfc 
Exact velccity magritude 



r.sqrt(|i-xO)-2.(y-yO)-2) 



Figure 4: 2D isentropic vortex. Initial data. Vortex velocity vectors (background velocity 
is substracted) and exact density, temperature, pressure and velocity magnitude as a 
function of r = \/ x'^ + y'^, {x' = x — xo,y' = y — yo) with (xq, yo) the center of the vortex. 



504 A vortex centered at Xq = (xo,yo) = (0,0) is added to the ambient gas at the initial 

505 time t = with the following conditions u = u^o + Su, v = Voo + Sv, and T* = + 6T* 



ou = —y — exp , ov = X — exp , 01 



(7-1)/? 2^ 



506 with r = \/ x'"^ + y'^, {x' = x — XQ,y' = y — yo) and vortex strength is given by /5 = 5.0. 

507 Consequently, the initial density is given by 



P = Pc 



T* \ -7-1 



1 



(7 - 
877r2 



exp (1 — r^) 



1 

7-1 



(51) 



508 and the pressure ig given by p = p"' . We assume periodic conditions and the exact solution 

509 at any time T > is the Sciine vortex but transla-ted. in the direction — (^005 -yoo). Note 

510 that Voo = (0, 0) generates a static vortex which is usually simpler to solve and can also 

511 be misleading. The exact density function for any point at time T is denoted by p'^^{x, T), 
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512 moreover in figure H] are plotted the exact solution for density, temperature, pressure and 

513 velocity as a function of r. A series of refined meshes (from 25 x 25 up to 400 x 400 cells) 

514 are successively used to compute the numerical solution where r = 10~^ is used. The 

515 meshes are made of 20^ points in velocity space with bounds [—15,15]^. The errors at 

516 time T = t"" for M spatial cells in one direction are given by 

517 The rates of convergence are computed as log(eM'/^Af )/log(M'/M) for two meshes with 

518 M' and M cells. In Table [1] are gathered the and L°° errors on the density variables 

519 and rates of convergence for FKS, the unlimited HOFKS and HOFKS at final time T = 1. 

520 As expected the FKS produces only first order accurate results in L} and L°° norms. 

521 Contrarily the unlimited HOFKS can reach a genuine higher order of accuracy in both 

522 norms; the high accuracy in L°° norm is due to the fact that the solution is smooth and no 

523 limiter is applied therefore extrema are only little diffused compared to limited schemes. 

524 Finally the (limited) HOFKS behaves, as expected, like a high order accurate scheme in 
norm and like a first order scheme in L°° norm. We also display in figure [5] the convergence 











Advected isentropic vortex 
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Unlim. 


HOFKS 




HOFKS 


Ax 


M 


err 


r oo 

Ij err 


err 


T OO 

Ij err 


err 


T OO 

Ij err 


1/25 




1.26E-02 




2.78E-01 




6.36E-03 




1.16E-01 




4.64E-03 




9.21E-02 




1/50 


502 


8.36E-03 


0.59 


2.07E-01 


0.43 


2.12E-03 


1.59 


3.77E-02 


1.62 


2.08E-03 


1.16 


5.60E-02 


0.72 


1/100 


100^ 


5.09E-03 


0.72 


1.22E-01 


0.76 


5.68E-04 


1.90 


l.lOE-02 


1.78 


6.40E-04 


1.70 


2.85E-02 


0.97 


1/200 


200^ 


2.86E-03 


0.82 


6.34E-02 


0.95 


1.49E-04 


1.93 


3.07E-03 


1.84 


1.64E-04 


1.97 


l.OOE-02 


1.51 


1/300 


300^ 


2.03E-03 


0.88 


4.27E-02 


1.01 


7.08E-05 


1.83 


1.55E-03 


1.69 


7.55E-05 


1.91 


6.56E-03 


1.04 


1/400 


400^ 


1.56E-04 


0.91 


3.17E-02 


1.03 


4.36E-05 


1.69 


9.85E-04 


1.57 


4.52E-05 


1.79 


4.93E-03 


0.99 


Expected order 


1 


1 


2 


2 


2 


1 



Table 1: L and L°° errors and convergence rates for the isentropic vortex problem with 
FKS, unlimited HOFKS and HOFKS schemes. 

525 

526 curves corresponding to the errors of Table [H 

527 Finally in figure [6] we present the density results obtained by the three schemes: FKS, 

528 unlimited HOFKS and HOFKS. The top panel presents the density on the 50 x 50 mesh for 

529 all schemes versus the exact solution. The bottom panels present the 50 x 50 and 100 x 100 

530 cell mesh results: the density as a function of r = \/ x'"^ + y'^, (x' = x — XQ{T),y' = 

531 y — yQ{T)), where xo{T) = xq + u T and yo{T) = y^ + v T are the exact coordinates of 

532 the vortex center at final time T, is plotted for all cells in the domain. Doing so we can 

533 measure the "convergence" of the results when a finer mesh is used, the excessive diffusion 

534 of the first order scheme and the tendency of undershooting of the unlimited scheme. 



_ maxjej |/3^^(xj,i") - 



(52) 
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535 5.3 2D Sod shock tube problem 

536 We consider now the 2D/2D Sod test on a square [0, 2] x [0, 2]. The velocity space is also 

537 a square with bounds —15 and 15, i.e. [—15, 15]^, discretized with A^^, = 20 points in each 

538 direction which gives 20^ points. The domain is divided into two parts, a disk centered 

539 at point (1, 1) of radius = 0.2 is filled with a gas with density pi = 1, mean velocity 

540 ul = Q and temperature Tl = 5, whereas the gas in the rest of the domain is initiated 

541 with PR = 0.125, UR = 0, Tr = 4. The final time is tgnai = 0.07. 

542 We report results for two different values of the Knudsen number r = 10"'^ and 

543 r = 10~^. For these regimes, we can appreciate the differences between the high order 

544 fast kinetic scheme and the fast kinetic scheme. For larger r, as expected, the solutions 

545 furnished by the two methods are very close and thus we do not show the figures. In the 

546 case of r = 10~^ we compare the results between the HOFKS and the FKS method with 

547 a first order and a second order MUSCL scheme for the compressible Euler equations. In 

548 the case of r = 10"'^ we compare the results between the HOFKS and the FKS method 

549 with a DSMC method for the Boltzmann-BGK equation. For this latter, we employed on 

550 average 100 particles per cell and we averaged the solution over 100 realizations. 

551 In Figure [71 we report the profiles fixing x = 1 for respectively the density, the mean 

552 velocity in the x-direction and in the y-direction and the temperature using a 200 x 200 

553 mesh for r = 10~^. In Figure [HI we report the same profiles for the same spatial position, 

554 i.e. X = 1, for the same macroscopic quantities but for a larger value of the Knudsen 

555 number: r = 10~^. In this latter case, for the velocity in the x-direction, we did not 

556 report the solution for the DSMC method because the number of particles employed does 

557 not permit to compute the solution with sufficient precision. Observe, in fact, that in this 

558 test case, the final value of the x- velocity is of the order of 10"'^, which is a value that due 

559 to the statistical fluctuations is very difficult to capture with DSMC methods. 

560 Moreover in Figure [7] bottom panels, we report some magnifications of the same profiles 

561 which permits to better appreciate the differences between the methods. We observe that, 

562 as in the ID case, the accuracy of the FKS method lies between the first and the second 
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Figure 6: 2D Vortex problem at t = 1. Density results obtained by the three schemes: 
FKS, unHmited HOFKS and HOFKS. Top: 3D view of the density on the 50 x 50 mesh 
versus the exact solution. Bottom: 50 x 50 (left panel) and 100 x 100 (right panel) cells 
meshes, density as a function of r = \l x'"^ + y'^, {x' = x — xo{T),y' = y — yo{T)), where 
xq(T) = xq + u T and yo{T) = yo + v T are the exact coordinates of the vortex center, is 
plotted for all cells. 



563 order spatial accuracy for small r. On the other hand, we cannot see the differences 

564 between the HOFKS scheme and the second order MUSCL scheme for r = 10~^ even 

565 with the magnifications. This is because the same MUSCL scheme is used for solving 

566 the compressible Euler equations and for constructing the HOFKS scheme. In the case 

567 r = 10~^, the HOFKS method, as in the ID case, gives sharper solutions with respect 

568 the FKS method. This is also the case for the DSMC method which exhibits a larger 

569 numerical diffusion with respect to our kinetic scheme. For larger r, the HOFKS, the 
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570 FKS and the DSMC almost collapse on the same line, if the number of particles is chosen 

571 sufficiently large, and thus we did not report any result. 

572 To conclude this test case study, we report, in tabled the CPU time T, as well as the 

573 CPU time per time cycle Tcyde, the CPU time per cycle per cell Tceii and the number of 

574 cycles needed to perform the computation for different meshes in space for respectively the 

575 FKS, the HOFKS and the DSMC method when the Knudsen number is r = 10~3. For the 

576 HOFKS and the FKS schemes the meshes are fixed in velocity and made of 20^ points. For 

577 the DSMC method we choose the same number of particles which has been used to produce 

578 the results shown in the previous Figures, i.e. 100 particles in average per cell and 100 

579 realizations. We computed the computational effort for these three schemes repeating the 

580 same test using different meshes in space ranging from = Ny = 25 to Nx = Ny = 200. 

581 The results of this analysis can be summarized as follows. The HOFKS method is around 

582 1.5 times more expensive than the FKS method for all the cases studied. This new scheme 

583 is however still very efficient: around 11.5 minutes to compute the solution on a 200 x 200 

584 spacial mesh for a kinetic equation on a mono-processor laptop. For comparison the DSMC 

585 method, which is also known to be a fast method, requires around 505 minutes. The ratio 
536 between the two methods is about 44. Moreover DSMC gives less accurate solutions as 

587 seen on the y-component of velocity which can not be accurately captured. 

588 5.4 2D Implosion problem 

589 Finally we consider a 2D/2D implosion problem on the square [0,2] x [0,2] discretized 

590 with 100^ points. As in the previous test, the velocity space is a square but with larger 

591 bounds, i.e. —20 and 20, which means [—20,20]^, discretized with A^^^, = 30 points in each 

592 direction which gives 40^ points. The domain is divided into two parts, a disk centered at 

593 point (1, 1) of radius = 0.2 is filled with a gas with density pL = 0.125, mean velocity 

594 = and temperature Tl = 4, whereas the gas in the rest of the domain is initiated 

595 with pn = 1, T/j = 4, velocity in the rr-direction = 1 for x E [0,1] and Ux = —1 for 

596 X € [1,2] while the velocity in the y-direction is initiated with Uy = 1 for y G [0, 1] and 

597 lij; = — 1 for y G [1, 2]. The final time is tfinai = 0.07. 

598 We report the results for the Knudsen number equal to r = 10"'^ comparing our scheme 

599 to the DSMC method for the Boltzmann-BGK equation. For this latter, we employed on 

600 average 200 particles per cell and the solution is averaged over 50 realizations. 

601 In FigureOwe report the isolines of density, x-velocity, y- velocity and temperature, for the 

602 HOFKS method on the left and the DSMC method on the right. We observe that the two 

603 methods furnish the same results except for the statistical noise of the DSMC method. On 

604 the other hand, the computational costs of the two approaches are still very different. We 

605 report, as for the 2D Sod test, in tabled the CPU time T, as well as the CPU time per 

606 time cycle Tcydc, the CPU time per cycle per cell Tceii and the number of cycles needed 

607 to perform the computation for different meshes in space. For the DSMC method we 

608 choose the same number of particles which has been used to produce the results shown in 

609 the figures, i.e. 200 particles in average per cell and 50 realizations. The results of this 

610 last test can be summarized as follows: The HOFKS method takes around 22 minutes for 

611 computing the solution on a 200 x 200 mesh for a kinetic equation on a mono-processor 
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Cell X # iVc 


# Deg. freedom Nt^t 


Cycle 


Time 


Time / cycle 


Time/cell 


X Ny 


N., xNyX Nl 


-'X:;yclc 


T(s) 


^cyclc (s) 


TccW (s) 


FKS scheme 


25 X 25 


25 X 25 X 20^ 


13 


- 1.5s 


0.12 


1.9 X 10""^ 


= 625 


= 250000 










50 X 50 


50 X 50 X 20^ 


25 


6s 


0.24 


1.04 X 10"^ 


= 2500 


= 10^ 










100 X 100 


100 X 100 X 20^ 


50 


50s 


1 


1.0 X 10""^ 


= 10000 


= 4 X 10^ 










200 X 200 


200 X 200 X 20^ 


100 


440s 


4.4 


1.1 X 10""^ 


= 40000 


= 16 X 10^ 










HOFKS scheme 


25 X 25 


25 X 25 X 20^ 


13 


~ 2s 


0.15 


~ 2.0 X 10""^ 


= 625 


= 250000 










50 X 50 


50 X 50 X 20^ 


25 


9s 


0.36 


1.44 X 10"^ 


= 2500 


= 10^ 










100 X 100 


100 X 100 X 20^ 


50 


77s 


1.54 


1.54 X 10"^ 


= 10000 


= 4 X 10^ 










200 X 200 


200 X 200 X 20^ 


100 


690s 


6.9 


1.72 X 10"^ 


= 40000 


= 16 X 10^ 










DSMC scheme 


25 X 25 


N(, X Naveragexcell 

25 X 25 X 100^ 


11 


73s 


6.63 


0.0106 


= 625 


= 6.25 X 10*^ 










50 X 50 


50 X 50 X 100^ X 50 X 50 


22 


540s 


24.54 


0.0098 


= 2500 


= 2.5 X 10^ 










100 X 100 


100 X 100 X 100^ 


45 


3700s 


82.22 


0.0082 


= 10000 


= 108 




~ 61mn 






200 X 200 


200 X 200 X 100^ 


90 


30300s 


336.66 


0.0084 


= 40000 


= 4 X 10^ 




~ 505mn 







Table 2: 2D Sod shock tube. Computational effort for the FKS, HOFKS and DSMC 
schemes for r = 10"'^. The tune per cycle is obtained by T^ycie = T/Ncyde and the time 
per cycle per ceh by Tceii = T/A'cycie/^c- 



612 laptop. The augmentation of the computational cost with respect to the 2D Sod test is 

613 essentially due to: first the augmentation of the mesh points in which the velocity space 

614 is discretized and second the reduction of the time step. The reduction of the time step is 

615 caused by the macroscopic solver of the compressible Euler equations which needs smaller 

616 time steps to ensure stability. For comparison the DSMC method, which still furnishes 

617 fluctuating and somewhat less accurate solutions, requires around 510 minutes. The ratio 

618 of CPU time is around 22 in favor of the HOFKS for a better overall accuracy. 
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Cell X # iVc 


# Deg. freedom 


Cycle 


Time 


Time/cycle 


Time/cell 


X Ny 




-^cycle 


T(s) 


^cycle (s) 


Tccll (s) 


HOFKS scheme 




N^xNyX Nl 










25 X 25 


25 X 25 X 30^ 


21 


~ 3s 


0.14 


~ 2.22 X 10^^ 


= 625 


= 562500 










50 X 50 


50 X 50 X 302 


43 


22s 


0.51 


2.04 X 10~^ 


= 2500 


= 2.25 X 10^ 










100 X 100 


100 X 100 X 30^ 


85 


180s 


2.12 


2.08 X 10-"^ 


= 10000 


= 9 X 10^ 










200 X 200 


200 X 200 X 30^ 


170 


1350s 


7.94 


1.98 X 10"^ 


= 40000 


= 36 X 10^ 




22.5mn 






DSMC scheme 


25 X 25 


Nx X Ny X N^veragexcell 

25 X 25 X 50 X 200 


9 


74s 


8.22 


0.0132 


= 625 


= 6.25 X 10^ 










50 X 50 


50 X 50 X 50 X 200 


20 


620s 


31 


0.0124 


= 2500 


= 2.5 X lO'^ 










100 X 100 


100 X 100 X 50 X 200 


40 


3842s 


96.05 


0.0096 


= 10000 


= 10^ 




~ 64mn 






200 X 200 


200 X 200 X 50 X 200 


81 


30600s 


377.77 


0.0094 


= 40000 


= 4 X 10^ 




~ 8.5h 







Table 3: 2D Implosion test. Computational effort for the HOFKS and the DSMC scheme. 
The time per cycle is obtained by Tcydc = T/A'^cycic and the time per cycle per cell by 

Tcell = T/A^cycle/^c- 
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density t=10 " 



velocity x=10 " 




Figure 7: 2D Sod test: solution at tgnai = 0.07 and x = I for r = 10~^. Top-Middle: Den- 
sity (top left), velocity in the x-direction (top right), velocity in the y-direction (middle left) 
and temperature (middle right). Bottom: magnification of the solution. HOFKS method 
continuous line, FKS method dash dotted line, first order and second order MUSCL meth- 
ods dotted lines. 
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0.2 0.4 




1.2 1.3 1.4 1.5 1.6 1.7 



temperature x=10 




1.4 1.5 1.6 1.7 




Figure 8: 2D Sod test: solution at tgnai = 0.07 and x = 1 for r = 10~^. Top- Middle: 
Density (top left), velocity in the x-direction (top right), velocity in the y-direction (middle 
left) and temperature (middle right). Bottom: magnification of the solution. HOFKS 
method continuous line, FKS method dash dotted line, DSMC dotted line. (DSMC results 
are not plotted for the y-component of velocity because the noise induced by the method 
produces oscillations the amplitude of which are far greater than the scale of the figure.) 
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Figure 9: 2D Implosion test at tfinai = 0.07 with r = IQ-^. HOFKS scheme (left), DSMC 
scheme (right). From top to bottom: density, temperature, x- velocity, y-velocity. 
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619 6 Conclusions 



620 In this work we have presented an high order in space extension of a new super efficient 

621 numerical method for solving kinetic equations. The method is based on a splitting be- 

622 tween the collision and the transport terms. The collision part is solved on a grid while 

623 the transport linear part is solved exactly by following the characteristics backward in 

624 time. The key point is that, conversely to semi-Lagrangian methods, we do not need to 

625 reconstruct the distribution function at each time step. This permits to tremendously 

626 reduce the computational cost with respect other existing methods for kinetic equations. 

627 In this paper, in order to solve the limitations in term of spatial accuracy close to the ther- 

628 modynamical equilibrium of the original Fast Kinetic Scheme, we coupled the solution of 

629 the FKS method with the solution of the compressible Euler equations. Then, we matched 

630 the moments obtained from the solution of the macroscopic equations with the moments 

631 obtained from the solution of the equilibrium part of the kinetic equation. Finally, we 

632 recovered the solution as a convex combination of the two contributions: the macroscopic 

633 and the microscopic parts. This improvement permits to preserve the desired accuracy in 

634 space for all the different regimes studied. 

635 The numerical results show that the HOFKS method performs as the FKS method 

636 for large values of the Knudsen number and as a high order shock capturing scheme for 

637 small Knudsen numbers. Moreover, the method requires a small computational effort. 

638 Numerical experiments have shown that the computational cost of the new method is 

639 around 1.5 times larger than the previous FKS method for a clear gain in accuracy when 

640 reached fluid regimes. Most importantly, we showed that this new class of fast kinetic 

641 schemes is more accurate and around 25 — 65 times faster than DSMC methods which are 

642 known to be efficient schemes. This important result opens the gate to extensive realistic 

643 numerical simulations of far from equilibrium physical models. 

644 In this work, we only focused on the spatial accuracy of the method and not on the 

645 time accuracy, we remind to future works for the development of schemes which are both 

646 accurate in time and in space. We also would like to extend the method to non uniform 

647 meshes and different discretizations of the velocity space. Finally, we want to apply this 

648 method to other kinetic equations as the Boltzmann or the Vlasov equation. 
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